function [tsdat,stepdat,erdat,fitdat]=procerode3(filein,filein2,crinfo)
%function [tsdat,stepdat,erdat,fitdat]=procerode3(filein,filein2,crinfo)
% mfile to process Gust microcosm erosion data and output in .mat and spreadsheet
% format
%
% filein = filename of datalogger file with turbidimeter data from erosion exp
%
% filein2 = filename of datafile with processed filter data from erosion exp
% for turbidimeter1
%
% temp= water temp for replacement water for turbidimeter1
%
% coren = core name
%
% crinfo = structured array with core specific information
%   crinfo.ntun = turbidimenter number (1 or 2)
%   crinfo.rpmn = motor number (1 or 2)
%   crinfo.subsamp = 'y' or 'n' indicating if bottle samples are subsamples
%       or effluent from entire step. If 'y' will read in file with "a"
%       calibration data
%
% nsteps = number of steps in experiment
%
% Aeff = effective area of sediment in core, [m^2]
%
% tstr = cell array of times for start of exp (tstr) and for
% end of ea step thereafter
%
% modified from procerode2 by PJD 26-Nov-2010
% - removed evals
% - simplify plotting
% - run one core at a time
% - add expfit subfunction for fitting, eliminating need for curve fitting
% toolbox
% - add code to read UGEMS files - must be switched in code to read
% original Gust files
% - add option of loading file with "a" and "q" rather than calculation for
% experiments where bottle samples are subsamples
% calc em by interating turbidity time series - allows subsampled and
% reqular experiments to be calculated in same way
% 
% 4-Apr-2011   bug fixed by PJD  - problem with indexing was occasionally
% causing inconsistent vector lengths.

%%
% read in data files
% [FILEID,YR,JDAY,HRMIN,SEC,NTU1,NTU2,RPM1,RPM2]=textread( filein,'%f%f%f%f%f%f%f%f%f%*[^\n]','delimiter',',');
% 
% [stepn,tau_nom,bottle_id,vol,sm,tss,ec,note]=textread(filein1,'%f%f%s%f%f%f%f%*s%s%*[^\r]','headerlines',1,'delimiter',',') %changed \n to \r
%  


% distribute variables from crinfo struct
titlestr = crinfo.titlestr;
temp = crinfo.temp;
coren = crinfo.coren;
nsteps = crinfo.nsteps;
Aeff = crinfo.Aeff;
tstr = crinfo.tstr;

% read in time series data file
fid1 = fopen(filein);

if(0)  % Gust system (e.g. Sanford's, Wiberg's) with CR10X datalogger
    C1 = textscan(fid1,'%f%f%f%f%f%f%f%f%f','delimiter',',');
    [FILEID,YR,JDAY,HRMIN,SEC,NTU1,NTU2,RPM1,RPM2] = C1{:};

    %find matlab time
    HR=floor(HRMIN/100);
    MM=mod(HRMIN,100);
    mtime=datenum(['31-dec- ' num2str(YR(1)-1)])+JDAY+HR/24+MM/(24*60)+SEC/(24*3600);
elseif(1)  % UGEMS system from GreenEyes
    
    C1 = textscan(fid1,'%s%f%f%f%f%f%f','delimiter',',');
    [dts,et,est,RPM1,NTU1,RPM2,NTU2] = C1{:};
    mtime = datenum(dts);
end
    

fclose(fid1);

% read in filter data file
if crinfo.subsamp == 'n'
fid2 = fopen(filein2);
C2 = textscan(fid2,'%f%f%s%f%f%s','headerlines',1,'delimiter',',');
[stepn,tau_nom,bottle_id,volml,smg,note] = C2{:};
vol = volml./1e6;
sm = smg./1000;

fclose(fid2);
elseif crinfo.subsamp == 'y'
    ssdat = load([coren '_a']);  % load file with a, q, etc calc from subsamples
    a = ssdat.a;             % calibration constant (g l-1 ntu-1)
    stepn = ssdat.stepn;     % step number
    sm = ssdat.sm;           % bottle sample sediment mass from filter (kg)
    q = ssdat.q;             % pumping rate (m3 s-1)
    vol = ssdat.vol;         % volume of bottle sample (m3)
    tau_nom = ssdat.tau_nom; % nominal applied stress (Pa)
    bottle_id = ssdat.bottle_id;  % bottle numbers of each step
    note = ssdat.note;            % note regarding bottle sample
end

tss = sm./vol;    % Total suspended solids (kg m-3) of effluent sample
etime=(mtime-datenum(tstr(1)))*24*3600;

% get time between individual datalogger samples (should be about 1 sec)
dt = [nan;diff(etime)];

% use ntun and rpmn inputs values to pick correct turbidity and rpm time
% series data
if crinfo.ntun == 1
    NTU = NTU1;
elseif crinfo.ntun == 2
    NTU = NTU2;
end

if crinfo.rpmn == 1
    RPM = RPM1;
elseif crinfo.rpmn ==2
    RPM = RPM2;
end

% trim data from before and after erosion expt.
idx.all = find(mtime >= datenum(tstr(1)) & mtime < datenum(tstr(end)) & ~isnan(NTU));
mtime=mtime(idx.all);
etime=etime(idx.all);
dt=dt(idx.all);
NTU = NTU(idx.all);
RPM = RPM(idx.all);

%loop to find indicies and times for ea step
for i1=1:nsteps
        
    idx.(['step' num2str(i1)]) = find(mtime >= datenum(tstr(i1)) & mtime < datenum(tstr(i1+1)));
    tstep(i1) = etime(idx.(['step' num2str(i1)])(length(idx.(['step' num2str(i1)]))))-etime(idx.(['step' num2str(i1)])(1));
end

TSS = [];

%loop to calc flow and TSS conversion for ea step
for i3=1:nsteps
      
    if crinfo.subsamp == 'n'
    q(i3) =vol(i3)./tstep(i3);
    a(i3) =sm(i3)./(q(i3).*sum(NTU(idx.(['step' num2str(i3)])).*dt(idx.(['step' num2str(i3)]))));
    end
    
    TSS = [TSS;NTU(idx.(['step' num2str(i3)])).* a(i3)];   
end

%find background TSS

% plot data from 1st step
fsz = 14;  % fontsize
set(0,'defaultaxesfontsize',fsz)
set(0,'defaulttextfontsize',fsz)

f_=figure
set(f_,'Units','Pixels','Position',[439.6 193.6 814.4 576.8]);
set(gcf,'name',['Background TSS Calc Core ' titlestr])

plot(etime(idx.step1),TSS(idx.step1),'b.','markersize',12)

xlim_(1) = min(etime(idx.step1));
xlim_(2) = max(etime(idx.step1));

% Nudge axis limits beyond data limits
if all(isfinite(xlim_))
   xlim_ = xlim_ + [-1 1] * 0.01 * diff(xlim_);
   xlim(xlim_)
end

xlabel('time [sec]');               % remove x label
ylabel('TSS [Kg ^. m^{-3}]');               % remove y label
set(gca,'yscale','log')

fsz = 14;  % fontsize
set(0,'defaultaxesfontsize',fsz)
set(0,'defaulttextfontsize',fsz)

f_=figure
set(f_,'Units','Pixels','Position',[439.6 193.6 814.4 576.8]);
set(gcf,'name',['Background TSS Calc Core ' titlestr])

plot(etime(idx.step1),TSS(idx.step1),'b.','markersize',12)

xlim_(1) = min(etime(idx.step1));
xlim_(2) = max(etime(idx.step1));

% Nudge axis limits beyond data limits
if all(isfinite(xlim_))
   xlim_ = xlim_ + [-1 1] * 0.01 * diff(xlim_);
   xlim(xlim_)
end

xlabel('time [sec]');               % remove x label
ylabel('TSS [Kg ^. m^{-3}]');               % remove y label
set(gca,'yscale','log')

% enter bounds of fit interval, fit, and plot
while cont ~= 'y'

    xmin=input('Enter lower bound of time [sec] for fit for Background TSS: ')
    xmin
    disp(' ')
    xmax=input('Enter upper bound of time [sec] for fit for Background TSS: ')
    xmax
    disp(' ')
    
    try
        [efit,hp1]=backgrnd_tss(etime(idx.step1),TSS(idx.step1),xmin,xmax);
        bckgrnd_TSS = efit(1).*exp(-efit(2).*4500);
        TSS_min = min(TSS);
        title(['Background TSS ' titlestr ' fit'])
        xlabel('time, [s]')
        ylabel('TSS, [kg ^. m^{-3}]')
        eqnstr=['TSS = ' num2str(efit(1)) ' ^. e^{-' num2str(efit(2)) ' * t}'];
        ylim=get(gca,'ylim');
        text(500,ylim(2)-0.2*(ylim(2)-ylim(1)),eqnstr)
        txtstr=['Background TSS = '  num2str(bckgrnd_TSS) ' [kg ^. m^{-3}]'];
        text(500,ylim(2)-0.4.*(ylim(2)-ylim(1)),txtstr)
        line([min(etime(idx.step1)) max(etime(idx.step1))],...
            [bckgrnd_TSS bckgrnd_TSS],'color','k');
    end

    disp(['Min measured TSS = ' num2str(min(TSS))])
    disp(' ')
    
    cont=input('Accept this value for background TSS? [y/n]: ','s')
    disp(' ')

% if fit not acceptable, enter new fit bounds or manually input background
% TSS value
    if cont == 'n'
        man=input('Enter Manually [y] ; or try fit method again [n]?: ','s')
        if man == 'y'
            bckgrnd_TSS =input(['Enter Background Concentration '...
                'for TSS [kg m^{-3}]: ']);
            %
            figure(f_)
            clf

            set(gcf,'name',['Background TSS Calc Core ' titlestr])
            plot(etime(idx.step1),TSS(idx.step1),'b.','markersize',12)
           xlim_(1) = min(etime(idx.step1));
            xlim_(2) = max(etime(idx.step1));

            % Nudge axis limits beyond data limits
            if all(isfinite(xlim_))
                xlim_ = xlim_ + [-1 1] * 0.01 * diff(xlim_);
                xlim(xlim_)
            end

            xlabel('time [sec]');               % remove x label
            ylabel('TSS [Kg ^. m^{-3}]');               % remove y label
            set(gca,'yscale','log')
            ylow = 10^floor(log10(min(TSS(idx.step1))));
            yhigh = 10^(floor(log10(min(TSS(idx.step1))))+3);
            set(gca,'ylim',[ylow yhigh])

            title(['Background TSS ' titlestr ' fit'])
            xlabel('time, [s]')
            ylabel('TSS, [kg ^. m^{-3}]')

            text(500,ylim(2)-0.2*(ylim(2)-ylim(1)),'Manually Selected')
            txtstr=['Background TSS = '  num2str(bckgrnd_TSS) ' [kg^. m^{-3}]']
            text(500,ylim(2)-0.4.*(ylim(2)-ylim(1)),txtstr)
            line([min(etime(idx.step1)) max(etime(idx.step1))],...
                [bckgrnd_TSS bckgrnd_TSS],'color','k')
        break
        end
    end

end


saveas(gcf,['TSS_backgrnfit_' coren],'fig')

%calc corrected TSS
TSS_corr=TSS - bckgrnd_TSS;

%calc Erosion rate

E1 = [];
    
for i4=1:nsteps
    
      % calculate mass eroded for each step by integrating TSS time series
      em(i4) = sum(TSS_corr(idx.(['step' num2str(i4)])).*q(i4).*dt(idx.(['step' num2str(i4)])))';
      
      % calculate erosion rate time series from TSS
      E1 = [E1;TSS_corr(idx.(['step' num2str(i4)])).*q(i4)./Aeff];
end

ec = cumsum(em);
%calc shear stress (tau) for ea. time point
u1_15=-1e-5.*RPM.^2+0.011.*RPM+0.0956; %[cm/s]
u1=u1_15-u1_15.*(temp-15).*0.03./5; %[cm/s]
tau=u1.^2./10; %[Pa]

%create structured data output

tsdat=struct('core',coren,'mtime',mtime,'etime',etime,'ntu', ...  
    NTU,'rpm',RPM,'dt',dt,'tss',TSS,'tss_corr',TSS_corr,'E',...
    E1,'tau',tau);

tsdat.fields = {'tsdat = time series data';
                'core = core name';
                'mtime = time vector in matlab datenum format';
                'etime = elapsed time (s)';
                'ntu = turbidity time series from Hach Turbidimeter (ntu)';
                'rpm = motor rpm time series';
                'dt = time between individual time series samples (s)';
                'tss = Total suspended solids time series (kg m-3)';
                'tss_corr = corrected TSS time series (kg m-3) = TSS - background';
                'E = erosion rate time series (kg m-2 s-1)';
                'tau = applied stress time series calculated from RPM and temp (Pa)'};
                
stepdat=struct('core',coren,'idx',idx,'tstep',tstep,'q',...      
    q,'a',a,'stepn',stepn,'tau_nom',tau_nom,'vol',...
    vol,'sm',sm,'em',em,'tss' ,tss,'ec',ec,'bottle_id',{bottle_id},'note',{note});

stepdat.fields = {'stepdat = data / info for each step';
                  'core = core name';
                  'idx = step specific indices for time series data';
                  'tstep = length of each step (s)';
                  'q = pumping rate (m3 s-1)'
                  'a = NTU to tss calibration (kg m-3 ntu-1)';
                  'stepn = step number';
                  'tau_nom = nominal applied stress (Pa)';
                  'vol = volume water collected (m3)';
                  'sm = mass of sediment collected (kg)';
                  'em = eroded mass (kg) = integrated TSS after removing background';
                  'tss = Total suspended solids (kg m-3) of effluent sample';
                  'ec = cumulative eroded mass (kg)';
                  'bottle_id = ID of effluent bottles collected for turbidity calibration';
                  'note = note regarding bottle sample or filtration, 0 = regular sample, 1 = subsample'};
 
plt_er1(tsdat,stepdat,coren,nsteps,titlestr);

[fitdat]= lambdas(tsdat,stepdat,coren,nsteps,titlestr);

plt_er(tsdat,stepdat,fitdat,coren,nsteps,titlestr)

pause

[erdat,hp1,hp2]=er_results(tsdat,stepdat,fitdat,titlestr,coren,nsteps,Aeff);

close('all')                 

save(coren,'crinfo','tsdat','stepdat','erdat','fitdat')
%--------------------------------------------------------------------------

function [efit,h_]=backgrnd_tss(X1_tss0,TSS1_1,xmin,xmax)
% [efit,h_]=backgrnd_tss(X1_tss0,TSS1_1,xmin,xmax)
%
% fit exponential to TSS time series from first step. used to determine
% background TSS in replacement water
%
% Inputs
% X1_tss0 = elapsed time
% TSS1_1 = Total suspended solids time series
% xmin,xmax = define start and end point of data to be used in fit
%
% Outputs
% efit = output fit values, efit function is in the form y=a*exp(-b*x)
% where efit(1) = a and efit(2) = b
% h_ = figure handle
% Set up figure to receive datasets and fits

f_ = clf;
figure(f_);
set(f_,'Units','Pixels','Position',[439.6 193.6 814.4 576.8]);
 
% --- Plot data originally in dataset "TSS1_1 vs. X1_tss0"
X1_tss0 = X1_tss0(:);
TSS1_1 = TSS1_1(:);
plot(X1_tss0,TSS1_1,'b.','MarkerSize',12)
hold on
xlim_(1) = min(X1_tss0);
xlim_(2) = max(X1_tss0);

% Nudge axis limits beyond data limits
if all(isfinite(xlim_))
   xlim_ = xlim_ + [-1 1] * 0.01 * diff(xlim_);
   xlim(xlim_)
end


% indices of data to include in fit
ftinds = ~(isnan(X1_tss0) | isnan(TSS1_1)) & (xmin<X1_tss0 & X1_tss0<xmax); 

% estimate starting points for fit using data
stp1 = TSS1_1(min(find(ftinds)));
stp2 = -(mean(log(TSS1_1(ftinds))) - log(stp1))./ floor(mean(find(ftinds)));

% fit exponential function to time vs TSS data
[efit,model] = expfit(X1_tss0(ftinds),TSS1_1(ftinds),[stp1 stp2]);

% Plot this fit
[sse, FittedCurve] = model(efit);
h_ = plot(X1_tss0(ftinds), FittedCurve, 'r','linewidth',1.5);
hold off;

xlabel('time [sec]');               % remove x label
ylabel('TSS [Kg ^. m^{-3}]');               % remove y label
set(gca,'yscale','log')
ylow=10^floor(log10(min(TSS1_1)));
yhigh=10^(floor(log10(min(TSS1_1)))+3);
set(gca,'ylim',[ylow yhigh])

function [estimates, model] = expfit(xdata, ydata,start_point)
    % exponential fit function using fminsearch
    % copied from Matlab help: Curve Fitting via Optimization
    % renamed, in help called fitcurvedemo, conveniently fittype was exactly
    % right for this purpose
    %
    % if fitting not producing good results, 'TolFun' and 'TolX' can be
    % altered. 
    %
    % added by Patrick Dickhudt 2-Dec-2010

    % Call fminsearch with starting points estimated from data values
    model = @expfun;
    options = optimset('TolFun',1E-6,'TolX',1E-6,'MaxFunEvals',1000,'MaxIter',1000);
    [estimates,fval,exitflag] = fminsearch(model,start_point,options);
    % expfun accepts curve parameters as inputs, and outputs sse,
    % the sum of squares error for A*exp(-lambda*xdata)-ydata,
    % and the FittedCurve. FMINSEARCH only needs sse, but we want
    % to plot the FittedCurve at the end.
    
function [sse, FittedCurve] = expfun(params)
        A = params(1);
        lambda = params(2);
        FittedCurve = A .* exp(-lambda * xdata);
        ErrorVector = FittedCurve - ydata;
        sse = sum(ErrorVector .^ 2);
        
end
end
end
end
%--------------------------------------------------------------------------